No evidence of DUI in the Mediterranean alien species Brachidontes pharaonis (P. Fisher, 1870) despite mitochondrial heteroplasmy

Two genetically different mitochondrial haplogroups of Brachidontes pharaonis (p-distance 6.8%) have been identified in the Mediterranean Sea. This hinted at a possible presence of doubly uniparental inheritance in this species. To ascertain this possibility, we sequenced two complete mitogenomes of Brachidontes pharaonis mussels and performed a qPCR analysis to measure the relative mitogenome copy numbers of both mtDNAs. Despite the presence of two very similar regions composed entirely of repetitive sequences in the two haplogroups, no recombination between mitogenomes was detected. In heteroplasmic individuals, both mitogenomes were present in the generative tissues of both sexes, which argues against the presence of doubly uniparental inheritance in this species.


Results
We were able to detect and sequence two slightly divergent mitogenomes in Brachidontes pharaonis samples (p-distance 0.071 for coding genes and 0.086 for the whole mtDNA). The preliminary results of the end-point PCR (Table 1) suggested not only the existence of heteroplasmy, but also showed that it was not linked with sex of individuals in this species. The results between PCR and qPCR differ because not every homoplasmic individual in PCR was checked with qPCR and some of the heteroplasmic individuals were resolved as homoplasmic after the qPCR. This was due to the very low estimated copy number of the second mitogenome, much lower than the copy number of a nuclear gene. Very high Cq values, close to the NTC control, were observed in these samples. Since such results could be either caused by non-specific signal (primer-dimer formation) or the truly very low gene copy number (much lower than that of atpα nuclear gene), these were classified as likely caused by sample contamination and not by true heteroplasmy. In the absence of any sex-bias, instead of naming those mitogenomes as M-type for male and F-type for female, we will refer to them neutrally as type A and type B. The mitogenomes (Fig. 1) are similar in length (type A: 20,066 bp, type B: 20,072 bp), code the typical set of 13 proteins, 2 rRNAs and 23 tRNAs, all on the same strand and in the same order. The striking feature of both mitogenomes is the presence of two long noncoding regions, rich in repetitive sequences, that divide the mitogenomes almost perfectly in half. Fragments of open reading frames in noncoding regions were identified as repeats of the 5′ end of atp6, suggesting the occurrence of multiple tandem duplication random loss events that encompass fragments of the coding sequences, throughout the evolution of these mitogenomes.
Given the overall relatively small divergence in mitochondrial sequences, the long repetitive sequences, and heteroplasmy, the possibility of widespread mitochondrial recombination at one of the long noncoding regions was considered. Such recombination would impair our ability to detect the type of mitogenome with just one small qPCR fragment. To address this issue, four pairs of specific qPCR primers were used to quantify the two parts of mitogenomes A and B separately (Supplementary data Table S1 and S2). Under the no recombination hypothesis the four quantities, one per part of the mitogenome, should match according to the mitogenome type, whereas if recombination is involved, the four estimates would be mismatched. Analysis of 30 B. pharaonis individuals revealed no noticeable difference between the halves of mitogenome copy number (Fig. 2). We have observed seven cases of heteroplasmy divided more or less evenly (no statistically significant differences for Fisher's exact test) between sexes. The heteroplasmic individuals were: 4 females, 3 males in qPCR (23%) and 6 females, 7 males in end-point PCR analysis (33%). Both A and B mitogenomes were detectable in males and females (Supplementary data Table S5-S10 and Fig. S1).
The newly reported mitogenomes of B. pharaonis form a monophyletic clade within the Brachidontes genus clade (Fig. 3). Brachidontes is a large genus composed of around 50 individual species (WoRMS database). However, due to the small number of available complete mitogenomes, the analysis of taxonomic affinities of the mitogenomes, including possible mitochondrial introgression events within this genus can not be conclusive. As a side note, phylogenetic analysis revealed the following inconsistency: genus Brachidontes seems to be paraphyletic, due to the presence of M. solisianus mitogenome within the clade. This hints that the taxonomy of M. solisianus may require reconsideration. Nevertheless, the phylogenetic distance between the two B. pharaonis mitogenomes is similar to the intraspecies distance between F-type mitogenomes form north and south populations of the Chilean mussels Perumytilus purpuratus which were separated from each other in Pleistocene [30][31][32] . It is also greater than the distance between several interspecies pairs of mitogenomes: Mytilus edulis and M. galloprovincialis, Bathmodiolus marisindicus and B. septemdierum, as well as all available mitogenomes form genus Gigantidas.

Discussion
We reported two mitogenomes (p-distance 8.6%) in B. pharaonis. These two mitogenomes are quite frequently present in the same individuals (Table 1). Is it possible that the observed heteroplasmy was not true but the result of a recent numt which is segregated in the population? We dismiss this possibility for the following reasons. www.nature.com/scientificreports/ Both mitogenomes were highly expressed at RNA level (31% of all NGS RNAseq sequencing reads mapped onto the A-type mtDNA and 4% onto the B-type mtDNA sequence), which is not a typical feature of a numt.  www.nature.com/scientificreports/ Moreover, numts are usually much smaller than the complete mitogenomes (but look here 35 ), and their sequence is degenerated 36 . Here both very large mitogenomes have all their genes intact and seemingly functional. Finally, the ratios of copy numbers of mitogenomes to the nuclear genome are not constant. A tight correlation between a numt and nuclear copy number is expected, which is not the case here. What can be the reason for the existence of two such haplotypes within this species? Since DUI was postulated for some Brachidontes species based solely on sequence divergence, it is reasonable to consider it here 3,24 . One may ask what sequence divergence threshold (if there is one) can be unequivocally associated with DUI? The situation present in the Baltic Sea Mytilus mussels 37,38 , where the divergent paternal mitogenome was replaced by the masculinized F-type mitogenome, is exceptional. Consequently, besides the differences in the noncoding control region, the M and F mitogenomes are identical in that case. If such extreme cases are not considered, the lowest divergence of the complete mitochondrial genomes observed and counted as DUI belongs to Arctica islandica 28,29 , with the M-F divergence (nucleotide p-distance) at 5.1% for coding genes and 6.9% for the whole mtDNA (Fig. S2). Even if a wider species range is considered, for which only cox1 16 gene fragments are known, the least divergent is also Arctica islandica: 2.5% at the protein level and 5.2% at the nucleotide level. So, can we count the divergence observed in B. variabilis 3 3.6% at the nucleotide level as a suggestion of DUI existence, with the lowest known divergence in the cox1 gene? Or is this just a high haplotype diversity within a species, not linked to sex-related heteroplasmy? This rhetorical question also applies to B. pharaonis.
In the context of DUI, divergence alone is just a secondary factor, depending only on the time of separation of the two lineages and their evolutionary rates. Nevertheless, for the emerging mitogenomes to play different physiological roles, there must be functional differences between their proteins, which is problematic for low genetic distances. However, even with a small number of nonsynonymous substitutions in the sequence, the encoded protein may become different enough (like in quickly evolving viruses or targeted enzyme engineering 39,40 ) to differ in activity. Unfortunately, tracing such substitutions in the context of DUI is difficult, but the possibility that a low-divergence mitogenome may still play distinct physiological role can't be discounted. Are there any other features that could help classify a mitogenome as involved in DUI? Perhaps we should consider additional open reading frames (ORFs) or gene extensions with no homology to the second mitogenome as a DUI marker? Is this a universal feature of DUI mitogenomes? Indeed, this seems to be true for all known pairs of DUI mitogenomes www.nature.com/scientificreports/ ( Table 2): each of them has some gender-specific structural features, but this is not the case for B. pharaonis.
Here, both mitogenomes are structurally identical. Another feature of mtDNA, which has been reported in the context of DUI, is mitochondrial recombination 65,66 . The existence of very large repetitive elements suggests, that at least intramolecular recombination would be present in B. pharaonis. Yet, we have not found any evidence for the exchange of large parts between A and B mitogenomes. They apparently maintain their integrity, despite the conditions, which should favor the exchange 67 . The qPCR method would not detect low frequency, short span recombination so this has to be considered cautiously. However, no detectable signature of recombination is present within the A and B mitogenomic sequences, therefore even if such recombination is present, it is limited to somatic cells.
Several studies postulated the existence of cryptic species within the Brachidontes genus (B. pharaonis/variabilis 3 and B. puniceus/exustus 24,68 ). These were usually argued by the distinct divergence (p-distance) at the nucleotide level (~ 7-20%) between mitochondrial sequences. However, at the protein level, these distances were much lower (p-distance ~ 0-1.5%). These must be interpreted cautiously because only short fragments of one mitochondrial gene are available (cox1: AY621835-AY621860; AY621862-AY621865; AY621909; AY621911; AY825105-AY825108; DQ836012; DQ836013; DQ836019-DQ836021) and the mitogenomic distances may be quite different (Fig. 4). However, all nonsynonymous substitutions in the mentioned cox1 fragments stay within the group of nonpolar, mostly aliphatic amino acids (Ile, Val, Leu, and in two cases also Met and Phe). This suggests no alternation in the overall protein structure/reactivity as is indicated by amino acid substitution matrices 69,70 derived from empirical data. The issue of cryptic species identification becomes even more complicated when specimens are erroneously assigned 68 to their respective species, as in the case of Brachidontes/   Orbigny, 1842). The current taxonomy classifies this species as Mytilaster, but the mitogenomic phylogenetic tree (Fig. 3) does not support this. Furthermore, the M. solisianus mitogenome deposited in GenBank has been wrongly annotated as Perna perna 71 . A taxonomic revision of Brachidontes seems warranted but should not be based solely on mitochondrial markers, which are known to evolve under strong selective constraints and are prone to introgression 72,73 . A good example of potential pitfalls associated with such a simplistic use of mitochondrial markers comes from the case of Baltic Mytilus trossulus, which in fact is a nuclear hybrid of M. edulis and M. trossulus that lost their native mitochondrial genome towards one from M. edulis 72,74,75 (also see here for M/F recombination 66,76 ). If only the mitogenomic protein level p-distances are considered, the cryptic species hypothesis is becoming less likely. However, the fact that heteroplasmic individuals were consistently observed suggests that anomalies in the mitochondrial inheritance may be involved.
In conclusion, B. pharaonis represents a species possessing two slightly different mitogenomes (p-distance 8.6%) in every possible combination between the sexes. Total homoplasmy for mitogenomes A or B, as well as heteroplasmy of both mitogenomes within a single individual are possible. This heteroplasmy is not correlated On a more practical level, we advise that the use of somatic tissues 77 , in phylogenetic studies on bivalves, usually practiced as a precaution to avoid amplification of the potential M mitogenome, may not always work as planned 78 . In the case of B. pharaonis such an approach would lead to seemingly random amplification of one of the two mitogenomes present in an individual.

Methods
Samples of Brachidontes pharaonis mussels were gathered in June 2014-2015 at the salt pan "infersa" of the Marsala lagoon (northwest of Sicily-Italy). Individuals were sectioned with a sterile scalpel blade, checked for male or female gametes under a light microscope, and stored frozen in − 80 °C until further use. DNA was extracted using the modified CTAB extraction method 79 Total RNA from three mantles of male individuals was pooled and sent to Macrogen Inc (Korea) for high throughput NGS sequencing (MiSeq Illumina, TruSeq NGS library 2 × 150 bp). Raw sequencing reads have been submitted to the SRA GenBank database under accession number SRR19141670. The acquired data were processed according to the Oyster river protocol 34,81 and assembled into the raw transcriptome. Long-range PCR primers for amplification of overlapping mitogenome fragments have been designed based on assembled contigs containing mitochondrial genes identified with Wise2 software 82 . PCRs were carried out in a volume of 20 µl, containing 25 ng of DNA, primers at 0.4 µM each, dNTP at 200 µM, 1.5 mM MgCl 2, and 0.4 U Phusion High-Fidelity polymerase (Thermo Scientific) suspended in GC buffer for difficult GC-rich templates. The PCR conditions were as follows: initial denaturation at 98 °C for 30 s; 35 cycles of denaturation at 98 °C for 10 s, annealing (Table S3 and S4) for 30 s and extension at 72 °C for 8 min. The final extension at 72 °C lasted 10 min. The amplified products were then assigned to their respective mitogenome (A and B) and sent for another NGS sequencing (MiSeq Illumina, TrueSeq NGS library, 2 × 300 bp). Complete mitochondrial DNA sequences have been recovered with NOVOplasty 83 software and validated by mapping NGS reads onto the assembled mitogenomes in CLC Genomics Workbench 9.5 (QIAGEN). The two acquired mitogenomes were then annotated with CRITICA 84 , Wise2 82 , GLIMMER 85 , ARWEN 86 , and nhmmer 87 . Mitogenome circular diagrams and compositional indices were drawn with MITOCONSTRICTOR 33,34 (https:// github. com/ aburz ynski/ mitoc onstr ictor). Annotated mitochondrial genomes for Brachidontes pharaonis were deposited in the GenBank under accession numbers ON464163 and ON464164.
Based on assembled transcriptome and two mitogenome sequences, a set of five qPCR primer pairs spanning: nuclear atpα gene and mitochondrial nad4, atp6 (from both mitogenomes; type A and type B) were designed in Primer3 88 . The specificity of each primer pair was verified, and there was no cross-amplification between A and B mitogenomic fragments. Quantitative PCR efficiency has been verified by running standard curves in nine repetitions with seven dilution points of samples. Reactions quantifying A and B mtDNA, as well as nuclear DNA, were performed on an ECO48 (Illumina/now PCRmax) Real-Time PCR System according to the qPCR kit (EurX) manufacturer instructions, in 10 µl reaction volume containing 1 × SG qPCR Master Mix, 2 µl of DNA (at ~ 15 ng/µl) and 0.5 µM of each primer. The thermal profile was as follows: initial denaturation at 95 °C for   (Table 3). All mitogenomes from the Mytilidae family available GenBank database (accessed May 2021) were used. The 12 mitochondrial protein coding genes (atp8 was omitted due to high divergence and in a few cases annotation problems 89 ) were used. Each gene was aligned separately, at protein level in MEGA7 90 software, with the ClustalW algorithm, to ensure that the codon structure of each gene is retained at the alignment level. Gap Extension and Gap Opening costs were set at 5, and the final alignments were visually inspected. No obvious alignment problems were encountered. Phylogenetic reconstruction was done using two approaches: Bayesian inference (BI) and Maximum Likelihool (ML). The same set of individual gene alignments (as separate data partitions) was used in both analyses.
BI was performed in BEAST2 91 . The parameters were as follows: nucleotide substitution model GTR + I + 4G, based on the results from the bModelTest package for BEAST2 software, relaxed log-normal clock, due to the varying evolution rates between F and M-type mtDNA and Yule prior to the common tree. The Markov chain was run in four replicates for 10 7 generations and sampled every 10,000th step. The convergence of samples was checked with Tracer 92 , the effective sample size for estimated parameters was greater than 200. ML was performed in IQ-TREE 93 , under the default parameters with ultrafast bootstrap approximation parameter set to 10,000 replicates. Substitution models for every partition were chosen with build-in ModelFinder tool. Model GTR + F + I + 4G was the best fitting for all protein genes except nd2, nd4, nd5 for which GTR + F + R5 was chosen.  www.nature.com/scientificreports/